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ABSTRACT 

The electron number density (Ne) distributions in solar chromosphere and corona are 
usually described with models of different nature: exponential for the former and inverse power 
law for the latter. Moreover, the model functions often have different dimensionality, e.g. the 
chromospheric distribution may depend solely on solar altitude, while the coronal number density 
may be a function of both altitude and latitude. For applications which need to consider both 
chromospheric and coronal models, the chromosphere-corona boundary, where these functions 
have different values as well as gradients, can lead to numerical problems. We encountered this 
problem in context of ray tracing through the corona at low radio frequencies, as a part of effort 
to prepare for the analysis of solar images from new generation radio arrays like the Murchison 
Widefield Array (MWA), Low Frequency Array (LOFAR) and Long Wavelength Array (LWA). 
We have developed a solution to this problem by using a patch function, a thin layer between the 
chromosphere and the corona which matches the values and gradients of the two regions at their 
respective interfaces. We describe the method we have developed for defining this patch function 
to seamlessly "stitch" chromospheric and coronal electron density distributions, and generalize 
the approach to work for any arbitrary distributions of different dimensionality. We show that 
the complexity of the patch function is independent of the stitched functions dimensionalities. It 
always has eight parameters (even four for univariate functions) and they may be found without 
linear system solution for every point. The developed method can potentially be useful for other 
applications. 



1. Introduction 

The visible surface of the sun, the photosphere, 
has the radius IRq ~ 6.955 x 10^ km. The pho- 
tosphere is surrounded by the solar atmosphere, 
which consists of two major layers. The lower 
layer, the chromosphere, is relatively thin. It ex- 
tends up to several thousand kilometers above the 
photosphere. The upper layer, the corona, is a 
gaseous envelope with the density rapidly decreas- 
ing with the radius. Depending upon the nature of 
the study, a working number for size of the corona 
can be determined based on the coronal height 
where the coronal plasma becomes too tenuous to 
be discernible for intended purposes. The size of 
the corona might thus be considered to range from 
tens up to hundreds of the solar radii. 

Electromagnetic rays at low frequencies in the 



heliosphere propagate in the regions with the 
plasma density lower than the critical density pcr, 
defined as 

= 

where e is the electron charge, rrie is the elec- 
tron mass, rup is the proton mass, and w is the 
radio wave frequency in rad s~^. At the crit- 
ical plasma density the dielectric permittivity e 
of the plasma and, hence, its index of refraction 
for a given frequency uj become zero. This im- 
plies that for most coronal models, radiation at 
frequencies < 150-180 MHz does not penetrate to 
the chromosphere and travels only in the corona. 
At higher frequencies, the rays with sharp angles 
of incidence can penetrate into the much denser 
and cooler chromosphere. Modeling of propaga- 
tion of radio frequency radiation through the so- 
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lar atmosphere hence requires knowledge of both 
chromospheric and coronal electron density (-/Ve) 
dist ributions. Our method o f choice for ray trac- 
ing ( Benkevitch et al.. 201Clt ). a fast, second order 



ability of high fidelity solar images from the new 
generation radio arr ays like the Murchison Wide- 



algorithm, requires smoothness in both and its 
gradient (VA^e) everywhere in the medium. The 
fact that chromospheric and coronal density dis- 
tributions come from independent models which 
are not consistent at the boundary between them 
leads to a discontinuity in and ViVg at this 
boundary. Though the transition between the 
chromosphere and the corona is intrinsically very 
rapid and involves rather large density gradients, 
the discontinuity in and ViVe is unphysical 
and leads to numerical instability in the ray trac- 
ing algorithm. It is therefore essential to make 
these distributions consistent at their interface and 
get them to smoothly morph from one to an- 
other while maintaining a smooth gradient. An 
attendant problem is the different dimensionality 
of the chromospheric and coronal density distri- 
butions: the former is one-dimensional, depend- 
ing only on the solar radius, while the latter can 
be two-dimensional, depending both on the radius 
and (co)latitude. 

In this work we offer our solution to this prob- 
lem. We use the concept of a thin patch region 
which encloses the chromosphere/corona interface, 
and which is consistent with the chromospheric 
and coronal models in either region. The patch 
is a function constructed to allow the two regions 
to smoothly morph into each other and accommo- 
dates the differences in their dimensionalities. In 
Section 2 we present the plasma density models 
and a mathematical formulation of the problem. 
In Section 3 we start with the simplest case of 
stitching two univariate functions. The method is 
developed further for the case of merging of a uni- 
variate and a multivariate functions in Section 4. 
The method appears to be universal enough to ex- 
tend it to the problem of merging two multivariate 
functions along one dimension. The appropriate 
theory is developed in Section 5. Numerical im- 
plementation details are discussed in Section 6, in 
particular it is shown that the suggested merging 
method is not computationally intensive. 

We came across this problem while building 
the formalism to implement ray tracing through 
the solar atmosphere at low radio frequencies. Our 
work is motivated by the expected near term avail- 



field Array (MWA) (jLonsdale et al.. II2009I) . The 



early results from the 32 element MWA engineer- 
ing prototype are very promising and currently 
represent the state of the art in hi gh fidelity so- 
lar im aging at these frequencies (joberoi et al..l 
20111 ). This work is of direct relevance for pur- 
suing solar science with other new generation 
arrays like the Low Frequency Array (LOFAR) 
( de Vos e t al.."2009') and Long Wavelength Array 
(LWA) (EUingson et al., 2009). We hope this for- 
mulation is general enough to be applicable to 
similar problems in other independent contexts. 

2. Formulation 

There are a number of models describing the 
electron number density 7Ve(cm~'^) in the chro- 
mosphere and the corona. For the chromo- 
sphere we use the mo del by Gillie and Menzel 
(ICillie fc Menzellll935l) : 



Ne{r) = 5.7 X ioiie-7-7xio-^(i?o(r-i)-500)^ 



(2) 



where r is the distance from the center of sun 
measured in the solar radii, i?©. This model is 
applicable for a chromosphere with a thickness of 
about 10,000 km, or O.O14378i?0. The simplest 
available coronal models are one dimensional in 
nature, depending o nly on the solar distance, e.g. 



the Newkirk model ( Newkirk 1 119611 ) 



iVeW =4.2 X 10'^+'^-32/'-, (3) 
and the Baumbach and Allen model ( Allen |[l947t l: 
Ne{r) = 10^(1. 55r-'5 + 2.99r"i6). (4) 

where r is the solar distance in units of Rq . 

However, observations show that the plasma 
density in the corona falls faster the solar distance 
in the polar regions than it does in equatorial re- 
gion, especially close to solar minima. This gave 
rise to development of more c omplex mod els. A 



popular one is by Kuniji Saito (jSaito Ill970l ). This 



coronal density model depends on two variables, 
the spherical coordinate r and 9, representing the 
radial distance from the center of the Sun and the 
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colatitude, respectively: 

Ne{r,e) = 3.09 X 10V"i6(l - 0.5 008 6*)+ 
1.56 X 10^r-^{l ~ 0.95 cos 6*)+ 
0.0251 X 10*r"2-^(l - Vcose). (5) 



i 
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Fig. 1. — Two functions /(r) and g{r,6), defined 
over different parts of space, are stitched together by a 
third, patch function p(r, 9) in the intermediate space 
segment between r = ri and r = r2. 

Consider a problem of stitching together the 
chromospheric and coronal models described in 
Eqs. ini and O respectively. Consider a rela- 
tively thin transition layer extending from ri to r2 . 
where ri is at the top of the chromosphere, and r2 
is at the bottom of the corona. For this specific 
problem the limits can be chosen as ri — 9, 000 km 
and r2 = 11,000 km above the photosphere. This 
transition layer encloses the chromosphere-corona 
boundary. Denote the iVe distribution functions 
in the chromosphere and the corona by f{r) and 
g{r, 9), respectively, as shown in Fig. [1] The patch 
function, p(r, 6), defined in the transition layer be- 
tween ri and r2, must match /(r) and g{r,9) at 
both its ends in both its value and derivatives. 
Finding an appropriate p(r,6) requires (1) choos- 
ing a family of functions for p(r, 6) and (2) finding 
its parameters to satisfy the boundary conditions 
for smoothness. 

Note that /(r) and g{r,9) have different num- 
ber of independent variables. Our approach 
can adequately accommodate this complication. 



though we start by considering a simpler case of 
stitching two one-dimensional distributions. 

3. Merging one- variable functions 

Consider the case where both chromospheric 
and coronal electron number density distributions 
are only dependent on one variable, the radial dis- 
tance r. As shown schematically in Fig.[Tl ri < r2, 
/(r) on the left is the Cillie and Menzel distribu- 
tion and g{r) on the right is either Newkirk 
or Baumbach distribution. We want the values of 
p{r) at ri and r2 to be equal to those of the dis- 
tributions on the left and right ends, respectively. 
Also, we want the derivative of p(r), denoted as 
p'{r), to take boundary values equal those of the 
derivatives of the distributions, f'(r) and g'{r) on 
the left and on the right, respectively. Putting the 
boundary conditions together in one system: 



P{ri) 


= f{ri) 


P{r2) 


= 9(^2) 


P'in) 


= .nri) 







(6) 



We have four equations, therefore p{r) must 
have exactly four unknown parameters to be fully 
determined through solving system ([5]). We also 
would like to have the simplest system to solve, 
i.e. a linear system. These considerations lead to 
the choice of a third order polynomial as p{r), 

p{r) = ao + air + a2r'^ + a^r^ . (7) 

Substitution of the polynomial ([7]) into the system 
yields the linear system 

' ao + airi + a2rl + asrf = /(ri) 
ao + air2 + a2rl + aar^ = gir2) 
ai + 2a2ri + Sa^rf = /'(ri) 
ai + 2a2r2 + 3a3r2 = .g'(?'2) 

II 11^ 
Its solution, the vector a = II ao; ai; 02; 03 1| , 

comprises the coefficients of polynomial ([7]) . Once 

these coefficients have been determined, the patch 
polynomial (O, by construction, has value and 
derivatives equal to those of the chromospheric 
distribution /(r) at its left end and the coronal dis- 
tribution g{r) at its right. The system ([S]) can be 
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Fig. 2. — Stitching of the two functions /(r) and g{r) 
by a patch function p(r) between r = ri and r = r2- 



also written compactly in the matrix-vector form, 
Pa = c, (9) 

II 1 1 T 

where c = ||/(ri); g{r2); f'{ri); g'{r2)\\ is 
the right-hand-side vector, and P is the system 
matrix. 
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(10) 



An example of smooth stitching of the chro- 
mospheric and coronal distributions, given by 
equations ([2]) and (jj]), is given in Fig. [2] 

4. Merging one- and multi-variable func- 
tions along one dimension 



number of constraints (or degrees of freedom) im- 
posed on the patch: two on the end point val- 
ues and four on the partial derivatives at the end 
points, or six on total. For a completely asym- 
metric coronal distribution g{r,9,ip), when it is a 
function of all the three variables, the radial dis- 
tance r, the colatitude 9, and the longitude ip, 
the number of conditions imposed on the patch 
p{r, 9, ip) becomes 2-1-3*2 = 8. These calculations 
suggest that the complexity of the patch function 
apparently grows with the dimensionality of the 
stitched functions. However, this is not the case. 
As we shall show here and in the next Section, the 
patch function can be build such that it always 
depends only on eight parameters independent of 
the dimensionality of the stitched functions. The 
method can be used for stitching functions of arbi- 
trary (and different) numbers of dimensions, along 
one of them. 

Note that the patch function must gradu- 
ally turn from multi-variable g{r,9,ip) at r2 into 
single- variable /(r) at ri. Consider a function 
g{r2,9,ip). It has its variable r "frozen" at r2, 
so that g{r2,9,ip) is now dependent only on the 
two remaining variables, 9 and ip. Form the prod- 
uct of the g{r2,9,Lp) and some polynomial b{r), 
which has value of at ri and 1 at r2. The func- 
tion b(r)g{r2, 9, ip) is not only equal to at ri, its 
partial derivatives with respect to 9 and tp also 
evaluate to at ri, i.e. at ri it effectively depends 
on r only. However, it cannot be linked to /(r) di- 
rectly because generally f{ri)y^O. To make such 
linkage possible, we add one more polynomial, 
a(r), so the patch function takes the form 

p{r, 9, ip) = a{r) + b{r)g{r2,9, ip). (11) 



We now generalize the solution to the case with 
a function of one variable at ri and a function of 
two or more variables at r2 . A real- world example 
here is the pro blem of merging the m odel of Gil- 
lie and Menzel (|Cillie fc Menzel Ill935l). Eg (H i n 
the chromosphere with the model of lSaito I ( 19701 ) . 
Eq. (O in the corona. The patch p{r, 9) is now 
a function of two variables. It must smoothly 
match both values and derivatives of /(r) at its 
left boundary and of g{r, 9) at the right bound- 
ary. Note that /(r) has only one derivative, while 
g{r, 9) has two partial derivatives. In order to 
guess the kind of a patch function we count the 



The polynomial a{r) must be equal f{r) at ri, 
and must be equal zero at r2, where b(r2)g{r2, 9, ip) 
equals g{r, 9, ip). Thus, the conditions imposed on 
the polynpmials are 



a{ri) = f{ri); 
bin) = 0; 



a(r2) = 0; 
6(r-2) = 1. 



(12) 



For the multi- variable functions the boundary con- 
ditions include their partial derivatives. Listing 
the boundary conditions as a system of equations 
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yields: 



for the derivatives of b{r), 



p{ri,e, ip) = a{ri) + b{ri)g{r2,0, if) = /(ri) 
Pr{ri,0,ip) = ar{ri) + br{ri)g{r2,0, ip) = fr{ri) 
Pe {ri,e,ip) ^ 6(ri )ge{r2,9,Lp) ^ Q 
{'ri,0,ip) ^ h{ri (rz , 0, = 
p{r2,0, (p) = a(r2) + b{r2)gir2,0, (p) 
= 9{r2,0,ip>) 
Pr{r2,0, v) = ar{r2) + br{r2)g{r2,0, ip) 

= 9r{r2,0,ip) 

Pe{r2,0, ip) = b{r2)ge{r2,0, ip) = ge{r2,d, f) 
^p^{r2,9,ip) = b{r2)g^{r2,0,(p) = g^{r2,0,ip) 

(13) 

Here we denote a partial derivative with respect to 
a variable by putting this variable at the subscript 
position after the function name, e.g. grif, (p) is 
the derivative of g{r, 9, ip) with respect to r. Now 
we simplify the system (|13p using conditions 
on the polynomials a{r) and b{r): 

p{ri,9,<p) = a(ri) = /(ri) 
Pr{ri,0,ip) = ar{ri) + briri)g{r2,0,ip) = fr{ri) 
pg{n,0, cp) = = 
p^{ri,9,^) = = 
p{r2,0, ip) = g{r2,9, ip) = g{r2,0, (p) 
Pr{r2,0, v) = ar{r2) + br{r2)g{r2,0, ip) 

= 9r{r2,0,'p) 
Pe{r2,0,ip) = g0{r2,d,p) = ge{r2,6,p) 
p^{r2,9,<p) = g^{r2,0,ip) = g^{r2,0,(p) 

(14) 

We can see that the first, the third, the fourth, the 
fifth, the seventh, and the eighth equations in (IT4)) 
are nothing but tautologies and can be removed 
from the system, which now takes the form 



ar(ri) + br{ri)g{r2,d,p) = fr{ri) 

= Pr{ri,e,ip) 
ar{r2) + br{r2)g{r2,0,ip) = gr{r2,0,ip) 

^Pr{r2,d,(p) 



(15) 



If we for simplicity fix the a{r) derivatives at zero, 
this set of equations provides boundary conditions 



(ariri) = 

ar{r2) = 

br{ri) 
br{r2) 



frin) 



9(^2, 0,ip) 

9r{r2,0,(p) 
9{r2,0,ip) ' 



(16) 



We remember that there are four constraints im- 
posed on each of the polynomials: two on their 
end-point values and two on their partial deriva- 
tives with respect to the stitching-dimension vari- 
able r. Therefore, both polynomials are of the 3'"'' 
order, have four unknown coefficients each, and 
require eight equations to be solved for. We have 
four equations in P^ : the other four are in (|12l) . 
Combining the both provides two linear systems 
that can be used to determine coefficients of the 
polynomials a{r) and b{r): 



a{ri) = f{ri) 

a{r2) = 

ar{ri) = 

. ar{r2) = 



(17) 



and 



b{r2) 

brin) 

br{r2) 



friri) 



9{r2,0,(p) 
9r{r2,0,'p) 



(18) 



9ir2,0,ip) ■ 

One can notice that both linear systems have a 
nice common property, they use the same sys- 
tem matrix P as used in the previous Section (see 
Eq. (HI through dTO])). The systems (HT]) and (HH]) 
can be rewritten in the vector-matrix form as fol- 
lows: 

Pa = Ca, (19) 
Pb = Cf,, (20) 



where a = 

ll&o; bi; 62, 
the polynomials 



|ao; ai; 

h 11^ 
03 



02; 03 1 1 and b = 
are vectors of coefficients for 



a{r) = ao + air + a2r'^ + asr^ . (21) 
b{r) = bn + bir + b2r^ + b3r\ (22) 
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Fig. 3. — Stitching of the f(r) function of a single 
variable and the function g{r,6) of two independent 
variables along the r dimension between r = ri and 
r — r2. The patch function p{r,6) has the form given 
by Eq. The g{r,0) function is plotted for three 

values of colatitude 9, 0°, 45°, and 90°. 



The right-hand-side vectors Cq and are deter- 
mined according to (|T7)) and ([18]): 



HI fin); 0; 0; ir, 



C6=|| 0; 1; 



g{r2,e,ip) ' 



gAr2,e,v) ||T 
g{r2,e,ip) I 



(23) 
(24) 



Substitution of a{r) and h{r) in pT|) gives a 
smooth patch function to link together the /() and 
g{) functions. 

In order to test the method we stitch two dis- 
tributions dependent on different number of vari- 
ables , those by Cillie a nd Menzel's ( Cillie fc Menzel 
19351) . Eq. © and bv lSaitol (|l97o|), Eq. The 
former depends only on r, while the latter depends 
on r and 9. The stitching occurs along the r co- 
ordinate. Fig. [3] shows three stitching cases for 
different colatitudes 9. This example illustrates 
the application of the general method described 
here. 

5. Merging two multivariable functions 
along one dimension 

It is possible to generalize the one-and-multi- 
variable merging method to stitch together two 
multi-variable functions along the axis of one of 
the variables they have in common. For example, 
we may need to smoothly merge two functions. 



f(r,9,ip) and g{r,9,Lp) along the r dimension us- 
ing the patch function p(r, 0, Lp) over the interval 
[ti , 7^2] . The patch function must have the proper- 
ties of /(r, 9, if) (value and derivatives) at the ri 
end, but on the way from ri to r2 these properties 
must smoothly morph into being indistinguishable 
from g{r, 9, ip). To build such a function, consider 
two functions formed from / and g by "freezing" 
the r variable: f{ri,9,if) and g{r2,9,ip). Follow- 
ing the method developed in the previous Section, 
we introduce two polynomials, a(r) and b{r). In 
the interval [ri,r2], a(r) should drop off from one 
to zero, while b{r) should increase from zero to 
one. Then the function: 

Pir,9,ip)^a{r)f{n,9,^) + bir)gir2,9,ip) (25) 

will behave as desired to serve as a valid patch 
function. The boundary conditions imposed on 
the values of the polynomials a(r) and b{r) are: 



(26) 



a{ri) = 1; a(r2) = 0; 
bin) - 0; 5(r2) = 1. 

As before, we build a system of equations compris- 
ing all the boundary conditions on p(r, 9, ip) at the 
end points ri and r2 similar to system 1)13^ : the 
values and derivatives of p{r, 9, ip) are equated to 
those of /(r, 0, ip) on the left and those oi g{r, 9, ip) 
on the right. After removing the tautologies only 
two equations remain: 

ar{ri)f{ri,9,ip) + br{ri)g{r2,9,ip) = fr{ri,d,ip) 

= Pr{ri,9,(p) 

ar{r2)f{ri,9, <p) + br{r2)g{r2,9, <p) = gr{r2,9, ip) 

= Pr{r2,9,(p) 
(27) 

This system of two equations has four unknowns, 
the derivatives of polynomials a{r) and b{r) at 
both ends of the interval [ri,r2]. For simplicity, 
we assign values of zero to ar{r2) and br{ri) to 
form the system 

fr{ri,9,ip) 



a-rin) 

ar{r2) 
br{ri) 

br{r2) 



fin, 




9r{r2, 



Kv) 



,v>) 



(28) 



9{'r2,9,(p) 

Combining conditions on the polynomials a(r) and 
b{r) from (1^ and (1^ we build two linear sys- 
tems: 
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and 



r a(^i) 

0(7-2) 

ar(ri) 

bin) 

br{ri) 
br{r2) 



1 


fr{ri,9,(p) 





1 


9r{r2, 



(29) 



(30) 



There are four constraints imposed on each of the 
polynomials, two on the end-point values and two 
on the partial derivatives with respect to r. Both 
polynomials are cubic, like (f2T|) and ([22|) . We al- 
ready have eight equations to solve for the eight 
coefficients of a(r) and b{r). Both system have 
the same system matrix P, given in (jlOp. and can 
be compactly written in the matrix-vector form 
as shown in ([T9l) and ([20l) . The right-hand-side 
vectors for systems ([29]) and (pO| are 



HI 1; 0; 



fririM.. 







Cb 



0; 1; 0; 



g('-2, 6,1,3) 



(31) 
(32) 



Thus a(r) and 6(r) are successfully found. The 
patch function (|25l) with these polynomials as co- 
efficients smoothly stitches together the two func- 
tions, f(r,0,(p) and g{r,9,ip), along the r dimen- 
sion on the interval [ri, r2]- Note that the method 
for stitching of two multivariable functions is not 
more complex than the method for stitching a 
single- variable function with a multivariable func- 
tion from the previous Section. 

So far we considered merging of the functions 
defined over real, three-dimensional space with the 
same set of independent variables. However, us- 
ing patch function in the form given in (j25p and 
the procedure described in this Section it is easy 
to show that the merged functions can have ar- 
bitrary numbers of variables, and even different 
numbers of variables. Having at least one variable 
in common is the only requirement. 

Consider a more general example. Suppose 
two functions defined in eight-dimensional space, 
/(x, y, z, g, r) and g{q, r, s, t, u) need to be smoothly 



merged along the q dimension over the interval 
[91,92]- The two functions have different sets of 
independent variables with non-empty intersec- 
tion having q as its element. The patch function 
must be dependent on the union of the both sets 
of variables. It is defined in the form 



p{x, y, z, q, r, t, u, s) = a{q)f{x, y, z, qi,r) 
+ b{r)g{q2,r,s,t,u). 



(33) 



Following the reasoning from the previous Sections 
we can show that the cubic polynomials a{q) and 
b{q) can be determined through solving the linear 
systems (fT9|) and (|20|) with the right hand side 
vectors 



Cb=|| 0; 1; 0- ^^^^^''•'^'''"^ 



g{q2 ,r,s,t,u) 



(34) 
(35) 



Discussion 



The method of stitching two functions along 
one dimension described here was developed in 
context of numerical simulations, which need to be 
optimized for computational efficiency and speed. 
If both chromospheric and coronal electron num- 
ber density distributions depend only on the dis- 
tance from the sun center r, the four coefficients 
of patch polynomial (O can be calculated at the 
initialisation step. During the simulation there is 
no need to solve system ([S]); only the polynomial 
values need to be calculated. For example, in the 
ray tracing application this only happens when the 
ray enters the transition layer between the chro- 
mosphere and the corona. Therefore, the com- 
putational overhead of the stitching is relatively 
small. 

In more complex cases of spherically non- 
symmetric coronal distributions like that by Saito 
([5]), the linear system ([20)1 must be solved at ev- 
ery step when the ray happens to travel inside the 
transition layer. This is required because the Cf, 
vector now depends on the 9 (and probably ip) 
variables, whose values may change from one step 
to the next. Here are a few tips to accelerate the 
computations. The fact that both linear systems 
have the same system matrix P ([T0| . is one of 
the benefits of the method. P only depends on 
the values of ri and r2 (i.e. boundaries of the 
transition layer), which are usually held constant 
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during the course of computations. This makes 
it possible to calculate P once and for all before 
the actual raytracing starts. To save even more 
time one can replace the linear system solving 
with matrix- vector multiplying. To do so, P~^, 
must be calculated during the initialization phase. 
The vector a of polynomial (f2T1) coefficients must 
also be pre-calculated as the product of the matrix 
P~^ and the vector Cq: 

a = P-^c,. (36) 

Then, in the process of ray tracing, only the coef- 
ficients of the polynomial b{r) (see (1^^ ) must be 
calculated as 

b = P-^Cb, (37) 

and followed by computing the patch function (fTTj) 
itself. Adopting this scheme makes computations 
significantly faster. Note, however, that in case of 
stitching two three-variable functions both oper- 
ations, ((36|) and (p6|) . must be repeated at every 
ray tracing step. 

In the previous Section it is shown that this ap- 
proach to stitching multi-variate functions is not 
limited to the three-dimensional cases. It is gen- 
eral enough to be applied to abstract functions of 
any number of independent variables, while retain- 
ing its simplicity. For any number of dimensions 
only eight parameters, the coefficients of the 3'"'^ 
order polynomials a{r) and b{r) need to be deter- 
mined. 

7. Conclusion 
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The method developed here is intended for 
smooth merging of the electron number density 
model distributions at the boundary between the 
solar corona and chromosphere. Used exten- 
sively in our coronal radio propagation studies, 
the method has shown excellent results. We have 
also shown that the method can easily be general- 
ized to the problems where two abstract functions 
of arbitrary variable numbers need to be smoothly 
stitched along one common dimension. 
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